Rana sierrae at Dusy Basin
Laetitia Wilkins
May 10 2019 in Munich
Background
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com. Sébastien Nusslé and I started this document to collaborate on our data analysis while he was in Switzerland and I was at Berkeley.
This is a project I have been working on with Kathleen Matthews, Zack Steel, Stephanie Carlson, and Sébastien Nusslé while I was in Stephanie Carlson’s lab https://nature.berkeley.edu/carlsonlab/. Kathleen had gone out every summer for many years to count yellow legged frogs at Dusy Basin in Kings Canyon National Park (Rana sierra, the Sierra Nevada yellow-legged frog). We compiled counts from 1997-2012. Volunteers went out at least three times a year (July, August, September) and counted the frogs at their various life-stages. We analyzed this dataset as a collaborative effort.
Acknowledgments
First, we would like to thank the summer field assistants (Adam Calo, Alex Figueroa, Ambrose Tuscano, Annie Shattuck, Carrie Sendak, Kyung Koh, Blake Meneken, Christopher Loomis, David Court, Ed Stauber, Elizabeth Gruenstein, Eva Patton, Igor Lacan, Jacob Martin, Jacquelyn Kracke, Jeff Schlueter, Jill Andersen, Joe Fontaine, Julia Anderson, Karen Pope, Krishna Feldman, Leah Culp, Linda Vance, Matthew Young, Michael Davis, Nick Earley, Rebecca Doubledee, Rick Ziegler, Robert Grasso, Steve Hasslinger, and Vincent King who did the surveys from 1997-2012. We would also like to thank Prof. Steve Beissinger (ORCID: 0000-0003-1323-2727) for his advice on the abundance analyses. Timothy Szeliga from NOAA’s Geo-Intelligence Division, William Brown from NOAA’s National Center for Environmental Information, and David Rizzardo from the California Department of Water Resources were especially helpful in finding publicly available precipitation data online. We are also thankful to Cassandra L. Ettinger (ORCID: 0000-0001-7334-403X), Sonia L. Ghose (ORCID: 0000-0001-5667-6876), and Marina E. De León (ORCID: 0000-0001-9973-1951) for their comments on earlier versions of this manuscript and the photograph of a Sierra Nevada yellow-legged frog in Fig. 1.
This script assumes you have a working directory in your home directory. Here, this directory is called “/DusyBasin”
All the input files for this script here are in “~/DusyBasin/Input_summaryStats/” or live: https://github.com/megaptera-helvetiae/DusyBasin/tree/master/Input_summaryStats
Set-up R and all the necessary libraries:
Take hash away if you still need to install these libraries!
# Take hash away if you still need to install these libraries!
#install.packages(corrplot)
#install.packages(scales)
#install.packages(lme4)
#install.packages(lmerTest)
#install.packages(unmarked)
library(corrplot)
library(scales)
library(lme4)
library(lmerTest)
library(unmarked)
p.trans <- function(p) {
res <- as.numeric(NA, length(p))
for (i in 1:length(p)) {
if (p[i] > 0.1) res[i] <- "p > 0.05 [NS]"
if (p[i] > 0.05 & p[i] <= 0.1) res[i] <- "p < 0.1 [.]"
if (p[i] > 0.01 & p[i] <= 0.05) res[i] <- "p < 0.05 [*]"
if (p[i] > 0.001 & p[i] <= 0.01) res[i] <- "p < 0.01 [**]"
if (p[i] <= 0.001) res[i] <- "p < 0.001 [***]"
}
return(res)}
# Data
load(file = "~/DusyBasin/Input_summaryStats/data.frogs.Rdata")First, we need a nice map for Figure 1!
Code for Figure 1:
# Since the end of 2018 we need an API key from Google to use their maps for ggmap. It is a hassle.
# You can learn more about it here:
# https://developers.google.com/maps/documentation/geocoding/get-api-key
# https://cloud.google.com/maps-platform/#get-started
# Once you signed up for it (with a credit card!!!):
# install dev version of ggmap
devtools::install_github("dkahle/ggmap")## Skipping install of 'ggmap' from a github remote, the SHA1 (c2c7f64d) has not changed since last install.
## Use `force = TRUE` to force installationlibrary(ggmap)## Loading required package: ggplot2## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.## Please cite ggmap if you use it! See citation("ggmap") for details.#> Loading required package: ggplot2
#> Google Maps API Terms of Service: http://developers.google.com/maps/terms.
#> Please cite ggmap if you use it: see citation("ggmap") for details.
# save api key
register_google(key = "AIzaSyCcOI7IYDFog73S-c8len4eEwmDd7FTWus")
# fill in your personal key!
bikemap2 <- get_map(location = c(-118.553, 37.103), maptype = "terrain", source = "google", zoom = 5)## Source : https://maps.googleapis.com/maps/api/staticmap?center=37.103,-118.553&zoom=5&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-c8len4eEwmDd7FTWusggmap(bikemap2)
bikemap3 <- get_map(location = c(-118.553, 37.103), maptype = "terrain", source = "google", zoom = 16)## Source : https://maps.googleapis.com/maps/api/staticmap?center=37.103,-118.553&zoom=16&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-c8len4eEwmDd7FTWusggmap(bikemap3)
# save it.Some more code to play with maps. It is not executed here. Change eval=TRUE, echo=TRUE and it will run.
Now let’s dive into the analysis!
Data exploration…
# Selection of meaningful lakes only (30 and above are rivers)
data <- subset(data, site %in% c(1:11, 13:29))
# Defining new categories
data$frogs <- data$subadults + data$adults
data$juveniles <- data$metamorphs + data$premetamorphs
data$tadpoles <- data$larvae + data$juveniles
plot(log(data$eggs+1) ~ log(data$larvae+1))
pred <- loess(log(data$eggs+1) ~ log(data$larvae+1))
points(pred$fitted[order(data$larvae)] ~ log(data$larvae+1)[order(data$larvae)], type = "l", lwd = 2, col = "red")
# New dataset: pooling seasons per year and sites as mean numbers
eg_mean_y <- tapply(data$eggs, list(data$year, data$site), mean, na.rm = TRUE)
la_mean_y <- tapply(data$larvae, list(data$year, data$site), mean, na.rm = TRUE)
ju_mean_y <- tapply(data$juveniles, list(data$year, data$site), mean, na.rm = TRUE)
su_mean_y <- tapply(data$subadults, list(data$year, data$site), mean, na.rm = TRUE)
ad_mean_y <- tapply(data$adults, list(data$year, data$site), mean, na.rm = TRUE)
fr_mean_y <- tapply(data$frogs, list(data$year, data$site), mean, na.rm = TRUE)
ta_mean_y <- tapply(data$tadpoles, list(data$year, data$site), mean, na.rm = TRUE)
ndat <- data.frame(year = as.numeric(rep(rownames(eg_mean_y), times = dim(eg_mean_y)[2])), site = as.numeric(rep(colnames(eg_mean_y), each = dim(eg_mean_y)[1])), eggs = as.numeric(eg_mean_y), larvae = as.numeric(la_mean_y), tadpoles = as.numeric(ta_mean_y), juveniles = as.numeric(ju_mean_y), subadults = as.numeric(su_mean_y), adults = as.numeric(ad_mean_y), frogs = as.numeric(fr_mean_y))
summary(ndat)## year site eggs larvae
## Min. :1997 Min. : 1.000 Min. : 0.0000 Min. : 0.00
## 1st Qu.:2000 1st Qu.: 4.000 1st Qu.: 0.0000 1st Qu.: 0.00
## Median :2005 Median : 7.500 Median : 0.0000 Median : 11.71
## Mean :2005 Mean : 9.214 Mean : 0.8973 Mean : 114.82
## 3rd Qu.:2009 3rd Qu.:11.000 3rd Qu.: 0.1250 3rd Qu.: 100.20
## Max. :2012 Max. :22.000 Max. :25.0000 Max. :2045.80
## NA's :19 NA's :19
## tadpoles juveniles subadults adults
## Min. : 0.00 Min. : 0.00 Min. : 0.000 Min. : 0.000
## 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.000 1st Qu.: 0.000
## Median : 11.71 Median : 0.00 Median : 0.800 Median : 2.250
## Mean : 137.95 Mean : 22.96 Mean : 6.748 Mean : 8.656
## 3rd Qu.: 131.00 3rd Qu.: 4.75 3rd Qu.: 4.589 3rd Qu.: 8.589
## Max. :2046.20 Max. :427.89 Max. :115.333 Max. :101.000
## NA's :19 NA's :19 NA's :19 NA's :19
## frogs
## Min. : 0.000
## 1st Qu.: 0.450
## Median : 3.333
## Mean : 15.446
## 3rd Qu.: 17.278
## Max. :170.000
## NA's :19Code for Figure 2
Let’s move on to proportion of occupied lakes:
# Year goodness of fit: proportion of lakes with eggs, larvae, juveniles, or adults
YGFeggs <- eg_mean_y
YGFeggs[eg_mean_y > 0] <- 1
YGFeggs <- apply(YGFeggs, 1, mean, na.rm = TRUE)
LGFeggs <- eg_mean_y
LGFeggs[eg_mean_y > 0] <- 1
LGFeggs <- apply(LGFeggs, 2, mean, na.rm = TRUE)
YGFlarvae <- la_mean_y
YGFlarvae[la_mean_y > 0] <- 1
YGFlarvae <- apply(YGFlarvae, 1, mean, na.rm = TRUE)
YGFjuveniles <- la_mean_y
YGFjuveniles[la_mean_y > 0] <- 1
YGFjuveniles <- apply(YGFjuveniles, 1, mean, na.rm = TRUE)
YGFadults <- ad_mean_y
YGFadults[ad_mean_y > 0] <- 1
YGFadults <- apply(YGFadults, 1, mean, na.rm = TRUE)
YGF <- rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults)
rownames(YGF) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
# First I have to add three more years with no counts to the plotting dataset:
more.years <- matrix(c(0,0,0,0,0,0,0,0,0,0,0,0), nrow=4, ncol=3, byrow=TRUE)
rownames(more.years) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
colnames(more.years) <- c("2013", "2014", "2015")
YGF <- as.matrix(YGF)
test <- cbind(YGF, more.years)
YGF <- test
par(mfrow = c(1,1), las = 2, mar = c(4,4,1,1))
barplot(YGF, beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(x = "topright", bty = "n"), col = terrain.colors(4), ylab = "Proportion of occupied lakes")
I also did an additional plot where I only included egg masses that were counted during the breeding season!
data.breeding <- subset(data, data$season=="b")
eg_mean_y_b <- tapply(data.breeding$eggs, list(data.breeding$year, data.breeding$site), mean, na.rm = TRUE)
YGFeggs <- eg_mean_y_b
YGFeggs[eg_mean_y_b > 0] <- 1
YGFeggs <- apply(YGFeggs, 1, mean, na.rm = TRUE)
YGFlarvae <- la_mean_y
YGFlarvae[la_mean_y > 0] <- 1
YGFlarvae <- apply(YGFlarvae, 1, mean, na.rm = TRUE)
YGFjuveniles <- la_mean_y
YGFjuveniles[la_mean_y > 0] <- 1
YGFjuveniles <- apply(YGFjuveniles, 1, mean, na.rm = TRUE)
YGFadults <- ad_mean_y
YGFadults[ad_mean_y > 0] <- 1
YGFadults <- apply(YGFadults, 1, mean, na.rm = TRUE)
YGF <- rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults)## Warning in rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults): number of
## columns of result is not a multiple of vector length (arg 1)rownames(YGF) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
# First I have to add three more years with no counts to the plotting dataset:
more.years <- matrix(c(0,0,0,0,0,0,0,0,0,0,0,0), nrow=4, ncol=3, byrow=TRUE)
rownames(more.years) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
colnames(more.years) <- c("2013", "2014", "2015")
YGF <- as.matrix(YGF)
test <- cbind(YGF, more.years)
YGF <- test
par(mfrow = c(1,1), las = 2, mar = c(4,4,1,1))
barplot(YGF, beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(x = "topright", bty = "n"), col = terrain.colors(4), ylab = "Proportion of occupied lakes")
And one with different colors:
YGFeggs <- eg_mean_y
YGFeggs[eg_mean_y > 0] <- 1
YGFeggs <- apply(YGFeggs, 1, mean, na.rm = TRUE)
LGFeggs <- eg_mean_y
LGFeggs[eg_mean_y > 0] <- 1
LGFeggs <- apply(LGFeggs, 2, mean, na.rm = TRUE)
YGFlarvae <- la_mean_y
YGFlarvae[la_mean_y > 0] <- 1
YGFlarvae <- apply(YGFlarvae, 1, mean, na.rm = TRUE)
YGFjuveniles <- la_mean_y
YGFjuveniles[la_mean_y > 0] <- 1
YGFjuveniles <- apply(YGFjuveniles, 1, mean, na.rm = TRUE)
YGFadults <- ad_mean_y
YGFadults[ad_mean_y > 0] <- 1
YGFadults <- apply(YGFadults, 1, mean, na.rm = TRUE)
YGF <- rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults)
rownames(YGF) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
# First I have to add three more years with no counts to the plotting dataset:
more.years <- matrix(c(0,0,0,0,0,0,0,0,0,0,0,0), nrow=4, ncol=3, byrow=TRUE)
rownames(more.years) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
colnames(more.years) <- c("2013", "2014", "2015")
YGF <- as.matrix(YGF)
test <- cbind(YGF, more.years)
YGF <- test
par(mfrow = c(1,1), las = 2, mar = c(4,4,1,1))
barplot(YGF, beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(x = "topright", bty = "n"), col = c("white", "green", "cornflowerblue", "coral1"), ylab = "Proportion of occupied lakes")
Now let’s do a combination of the two plots where I take the sums for egg mass but otherwise the data from the first figure with averages over the three sampling efforts each summer:
# for eggs take only breeding season:
data.breeding <- subset(data, data$season=="b")
eg_mean_y_b <- tapply(data.breeding$eggs, list(data.breeding$year, data.breeding$site), mean, na.rm = TRUE)
YGFeggs <- eg_mean_y_b
YGFeggs[eg_mean_y_b > 0] <- 1
YGFeggs <- apply(YGFeggs, 1, mean, na.rm = TRUE)
YGFlarvae <- la_mean_y
YGFlarvae[la_mean_y > 0] <- 1
YGFlarvae <- apply(YGFlarvae, 1, mean, na.rm = TRUE)
YGFjuveniles <- la_mean_y
YGFjuveniles[la_mean_y > 0] <- 1
YGFjuveniles <- apply(YGFjuveniles, 1, mean, na.rm = TRUE)
YGFadults <- ad_mean_y
YGFadults[ad_mean_y > 0] <- 1
YGFadults <- apply(YGFadults, 1, mean, na.rm = TRUE)
YGF <- rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults)## Warning in rbind(YGFeggs, YGFlarvae, YGFjuveniles, YGFadults): number of
## columns of result is not a multiple of vector length (arg 1)rownames(YGF) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
# First I have to add three more years with no counts to the plotting dataset:
more.years <- matrix(c(0,0,0,0,0,0,0,0,0,0,0,0), nrow=4, ncol=3, byrow=TRUE)
rownames(more.years) <- c("Eggs", "Tadpoles", "Subadults", "Adults")
colnames(more.years) <- c("2013", "2014", "2015")
YGF <- as.matrix(YGF)
test <- cbind(YGF, more.years)
YGF <- test
YGF[1,2]<- 0.7
YGF[1,13]<- 0
YGF[1,14]<- 0
YGF[1,15]<- 0
par(mfrow = c(1,1), las = 2, mar = c(4,4,1,1))
barplot(YGF, beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(x = "topright", bty = "n"), col = terrain.colors(4), ylab = "Proportion of occupied lakes")
Figures 3 and 4 were done using code in another script by Zack Steel!!!
Snow Data
snow2 <- subset(season_length, season == "winter")[,-4]
names(snow2) <- c("w.start", "w.end", "year", "w.dur")
snow2$year <- 1989:2016
snow3 <- subset(season_length, season == "summer")[,-4]; names(snow3) <- c("s.start", "s.end", "year", "s.dur")
names(snow3) <- c("s.start", "s.end", "year", "s.dur")
snow3$year <- 1989:2017
snow2 <- merge(snow2, snow3, all = TRUE)
snow2$snow.max <- tapply(dsnow$watcont, dsnow$year, max, na.rm = TRUE)[-1]
par(mfrow = c(2,2), las = 1, mar = c(3,4,1,1))
plot(watcont ~ date, dsnow, type = "l", ylab = "Water Content Equivalent")
plot(tapply(dsnow$watcont, dsnow$year, max, na.rm = TRUE) ~ c(1988:2017), type = "b", pch = 17, ylab = "Maximal snow")
plot(duration ~ year, subset(season_length, season == "winter"), type = "b", pch = 16, ylab = "Winter duration")
abline(lm(duration ~ year, subset(season_length, season == "winter")))
plot(duration ~ year, subset(season_length, season == "summer"), type = "b", pch = 17, ylab = "Summer duration")
abline(lm(duration ~ year, subset(season_length, season == "summer")))
Code for Figure 5
Snow plots
par(las = 1, mar = c(5,4,1,1), mfrow = c(2,1))
plot(SnHe_cm ~ Year, snow, type = "b", pch = 16)
boxplot(SnHe_cm ~ Year, snow, type = "b", pch = 16)
# Snow resampling (1 value per year)
boot = 100
years <- unique(snow$Year)
bootsnow <- matrix(NA, boot, length(years))
delta_snow <- NA
t0 <- Sys.time()
for (i in 1:boot) {
for (j in 1:length(years)) {
interm <- subset(snow, Year == years[j])$SnHe_cm
if (length(interm) > 1) bootsnow[i,j] <- sample(x = as.numeric(subset(snow, Year == years[j])$SnHe_cm), size = 1)
else bootsnow[i,j] <- interm
}
res <- lm(bootsnow[i,] ~ years)
delta_snow[i] <- coef(res)[2]
}
Sys.time()-t0## Time difference of 1.197238 secshist(delta_snow)
mean(delta_snow); sd(delta_snow)## [1] -0.6067818## [1] 0.1246391t.test(delta_snow)##
## One Sample t-test
##
## data: delta_snow
## t = -48.683, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.6315129 -0.5820507
## sample estimates:
## mean of x
## -0.6067818# Snow averaging (1 value per year)
snow.year <- tapply(snow$SnHe_cm, snow$Year, mean)
snow.year2 <- tapply(snow$SnHe_cm, snow$year2, mean)
snow.min <- tapply(snow$SnHe_cm, snow$Year, min)
snow.min2 <- tapply(snow$SnHe_cm, snow$year2, min)
snow.max <- tapply(snow$SnHe_cm, snow$Year, max)
snow.max2 <- tapply(snow$SnHe_cm, snow$year2, max)
par(las = 1, mar = c(4,4,1,1))
plot(snow.year ~ as.numeric(names(snow.year)), pch = 16, type = "n", ylim = c(0, 500), ylab = "Snow height [cm]", xlab = "Time [year]")
abline(v = c(194:202)*10, lty = 2)
polygon(c(years, rev(years)), c(snow.min, rev(snow.max)), col = "grey")
points(snow.year ~ as.numeric(names(snow.year)), pch = 16, type = "b")
abline(h = 100)
arrows(1961.5, 80, 1974.5, 80, lwd = 2, code = 3)
arrows(1977.5, 80, 1986.5, 80, lwd = 2, code = 3)
arrows(1992.5, 80, 2000.5, 80, lwd = 2, code = 3)
arrows(2001.5, 80, 2006.5, 80, lwd = 2, code = 3)
arrows(2007.5, 80, 2011.5, 80, lwd = 2, code = 3)
text(c(1967.5, 1981.5, 1996, 2004, 2010), y = 70, c(13,9,8,5, 4))
snowheight <- data.frame(years, snow.year, snow.min, snow.max)
res <- lm(snow.year ~ years, subset(snowheight, years %in% 1958:2015));summary(res)##
## Call:
## lm(formula = snow.year ~ years, data = subset(snowheight, years %in%
## 1958:2015))
##
## Residuals:
## Min 1Q Median 3Q Max
## -139.26 -62.46 -15.95 50.85 198.11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1655.6963 1276.5364 1.297 0.200
## years -0.7394 0.6426 -1.151 0.255
##
## Residual standard error: 81.92 on 56 degrees of freedom
## Multiple R-squared: 0.0231, Adjusted R-squared: 0.005654
## F-statistic: 1.324 on 1 and 56 DF, p-value: 0.2547res <- lm(snow.min ~ years, subset(snowheight, years %in% 1958:2015));summary(res)##
## Call:
## lm(formula = snow.min ~ years, data = subset(snowheight, years %in%
## 1958:2015))
##
## Residuals:
## Min 1Q Median 3Q Max
## -114.814 -56.987 -4.765 48.538 203.331
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1153.7636 1100.5153 1.048 0.299
## years -0.5074 0.5540 -0.916 0.364
##
## Residual standard error: 70.63 on 56 degrees of freedom
## Multiple R-squared: 0.01476, Adjusted R-squared: -0.002834
## F-statistic: 0.8389 on 1 and 56 DF, p-value: 0.3636res <- lm(snow.max ~ years, subset(snowheight, years %in% 1958:2015));summary(res)##
## Call:
## lm(formula = snow.max ~ years, data = subset(snowheight, years %in%
## 1958:2015))
##
## Residuals:
## Min 1Q Median 3Q Max
## -157.59 -75.55 -21.54 81.39 239.97
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2002.1621 1506.2239 1.329 0.189
## years -0.8956 0.7582 -1.181 0.242
##
## Residual standard error: 96.67 on 56 degrees of freedom
## Multiple R-squared: 0.02431, Adjusted R-squared: 0.006889
## F-statistic: 1.395 on 1 and 56 DF, p-value: 0.2425res <- lm(snow.year2 ~ as.numeric(names(snow.year2)));summary(res)##
## Call:
## lm(formula = snow.year2 ~ as.numeric(names(snow.year2)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -73.070 -22.493 0.649 22.771 55.006
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1644.0277 757.5098 2.170 0.0464 *
## as.numeric(names(snow.year2)) -0.7301 0.3835 -1.904 0.0763 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 38.73 on 15 degrees of freedom
## Multiple R-squared: 0.1946, Adjusted R-squared: 0.1409
## F-statistic: 3.624 on 1 and 15 DF, p-value: 0.07632res <- lm(snow.min2 ~ as.numeric(names(snow.year2)));summary(res)##
## Call:
## lm(formula = snow.min2 ~ as.numeric(names(snow.year2)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -68.524 -15.743 2.045 19.927 106.525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3829.7157 803.7712 4.765 0.000251 ***
## as.numeric(names(snow.year2)) -1.8863 0.4069 -4.635 0.000324 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 41.1 on 15 degrees of freedom
## Multiple R-squared: 0.5889, Adjusted R-squared: 0.5615
## F-statistic: 21.49 on 1 and 15 DF, p-value: 0.0003236res <- lm(snow.max2 ~ as.numeric(names(snow.year2)));summary(res)##
## Call:
## lm(formula = snow.max2 ~ as.numeric(names(snow.year2)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -157.22 -32.08 -10.32 39.00 147.94
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1585.8358 1344.8631 -1.179 0.257
## as.numeric(names(snow.year2)) 0.9728 0.6809 1.429 0.174
##
## Residual standard error: 68.77 on 15 degrees of freedom
## Multiple R-squared: 0.1198, Adjusted R-squared: 0.06111
## F-statistic: 2.041 on 1 and 15 DF, p-value: 0.1736plot(snow.year2 ~ as.numeric(names(snow.year2)), pch = 16, type = "b", ylim = c(0, 500))
points(snow.min2 ~ as.numeric(names(snow.year2)), type = "l")
points(snow.max2 ~ as.numeric(names(snow.year2)), type = "l")
#-----------------------------------
# Graphs
#-----------------------------------
boxplot(SnHe_cm ~ Year, subset(snow, Year > 1900))
abline(h = 100, col = "red")
#quartz("", 10, 7)
par(mfrow = c(1,1), las = 1, mar = c(2,4,0,1), oma = c(1,1,1,1))
plot(SnHe_cm ~ Time, snow, type = "n", ylab = "Snow cover [cm]", pch = 16)
abline(h = c(1:10)*50, lty = 2, col = "grey")
abline(v = c(194:202)*10, lty = 2, col = "grey")
points(SnHe_cm ~ Time, snow, type = "p", ylab = "Snow cover [cm]", pch = 16)
points(lowess(snow$SnHe_cm)$y ~ snow$Time, col = "red", lwd = 2, type = "l")
res <- lm(SnHe_cm ~ Time, snow); summary(res)##
## Call:
## lm(formula = SnHe_cm ~ Time, data = snow)
##
## Residuals:
## Min 1Q Median 3Q Max
## -164.12 -65.24 -17.59 60.73 273.98
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1330.3285 637.8121 2.086 0.0382 *
## Time -0.5717 0.3216 -1.778 0.0769 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 92.19 on 209 degrees of freedom
## Multiple R-squared: 0.0149, Adjusted R-squared: 0.01018
## F-statistic: 3.161 on 1 and 209 DF, p-value: 0.07689abline(res)
legend("topleft", p.trans(coef(summary(res))[2,4]), bty = "n")
par(mfrow = c(3,1), las = 1, mar = c(2,4,0,1), oma = c(1,1,1,1))
plot(snow.year ~ as.numeric(names(snow.year)), type = "b", ylab = "Mean snow cover")
points(lowess(snow.year)$y ~ as.numeric(names(snow.year)), col = "red", lwd = 2, type = "l")
res <- lm(snow.year ~ as.numeric(names(snow.year))); summary(res)##
## Call:
## lm(formula = snow.year ~ as.numeric(names(snow.year)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -140.93 -58.35 -11.23 43.64 197.29
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1359.0107 722.4623 1.881 0.0636 .
## as.numeric(names(snow.year)) -0.5885 0.3654 -1.611 0.1113
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 76.89 on 79 degrees of freedom
## Multiple R-squared: 0.03179, Adjusted R-squared: 0.01953
## F-statistic: 2.594 on 1 and 79 DF, p-value: 0.1113abline(res)
legend("topleft", p.trans(coef(summary(res))[2,4]), bty = "n")
plot(snow.min ~ as.numeric(names(snow.year)), type = "b", ylab = "Min snow cover")
points(lowess(snow.min)$y ~ as.numeric(names(snow.year)), col = "red", lwd = 2, type = "l")
res <- lm(snow.min ~ as.numeric(names(snow.year))); summary(res)##
## Call:
## lm(formula = snow.min ~ as.numeric(names(snow.year)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -129.419 -51.061 -3.118 44.789 192.067
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2269.2859 656.1875 3.458 0.000879 ***
## as.numeric(names(snow.year)) -1.0643 0.3319 -3.207 0.001939 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 69.84 on 79 degrees of freedom
## Multiple R-squared: 0.1152, Adjusted R-squared: 0.104
## F-statistic: 10.28 on 1 and 79 DF, p-value: 0.001939abline(res)
legend("topleft", p.trans(coef(summary(res))[2,4]), bty = "n")
plot(snow.max ~ as.numeric(names(snow.year)), type = "b", ylab = "Max snow cover")
points(lowess(snow.max)$y ~ as.numeric(names(snow.year)), col = "red", lwd = 2, type = "l")
res <- lm(snow.max ~ as.numeric(names(snow.year))); summary(res)##
## Call:
## lm(formula = snow.max ~ as.numeric(names(snow.year)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -146.41 -63.85 -22.69 58.52 255.04
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 482.8551 843.0128 0.573 0.568
## as.numeric(names(snow.year)) -0.1317 0.4264 -0.309 0.758
##
## Residual standard error: 89.72 on 79 degrees of freedom
## Multiple R-squared: 0.001206, Adjusted R-squared: -0.01144
## F-statistic: 0.09539 on 1 and 79 DF, p-value: 0.7582abline(res)
legend("topleft", p.trans(coef(summary(res))[2,4]), bty = "n")
Now the code for Table 1
setwd("/")
datab<-read.delim("sc.orig.txt", header=T)
str(datab)
datab$lake<-as.factor(datab$lake)
library(lme4)
# Models with season length
res<-glmer(dry ~ surf_mean + sum_len + (1|lake), data=datab, family=binomial)
# Need to scale mean surface area measures because they are not normally distributed!
resa<-glmer(dry ~ surf_lm + sum_len + (1|lake), data=datab, family=binomial)
summary(res)
summary(resa)
# Which model is better, surface area or depth?
resb<-glmer(dry ~ depth + sum_len + (1|lake), data=datab, family=binomial)
summary(resb)
anova(resa, resb)
# Depth models for table:
resb.2<-glmer(dry ~ depth + (1|lake), data=datab, family=binomial)
resb.3<-glmer(dry ~ sum_len + (1|lake), data=datab, family=binomial)
anova(resb, resb.2)
anova(resb, resb.3)
# Depth - snow pack models for table:
resf<-glmer(dry ~ snow_max + depth + (1|lake), data=datab, family=binomial)
resf.1<-glmer(dry ~ depth + (1|lake), data=datab, family=binomial)
resf.2<-glmer(dry ~ snow_max + (1|lake), data=datab, family=binomial)
anova(resf, resf.1)
anova(resf, resf.2)
# Surface area models for table:
resa<-glmer(dry ~ surf_lm + sum_len + (1|lake), data=datab, family=binomial)
resa.1<-glmer(dry ~ sum_len + (1|lake), data=datab, family=binomial)
resa.2<-glmer(dry ~ surf_lm + (1|lake), data=datab, family=binomial)
anova(resa, resa.1)
anova(resa, resa.2)
# Surface area - snow pack models for table:
resg<-glmer(dry ~ snow_max + surf_lm + (1|lake), data=datab, family=binomial)
resg.1<-glmer(dry ~ surf_lm + (1|lake), data=datab, family=binomial)
resg.2<-glmer(dry ~ snow_max + (1|lake), data=datab, family=binomial)
anova(resg, resg.1)
anova(resg, resg.2)Now my favorite plot; Fig. S1 in the manuscript
The summary plot using bubbles to show abundance of different life stages each each. Sites are placed according to their geographic location in the plot.
d <- as.dist(dist)
mds.coor <- cmdscale(d)
mds.coor <- mds.coor[c(1,2,4,5,6,7,9,13),]
par(mar = c(0,0,0,0), mfrow = c(3,5))
for (i in 1:dim(ad_mean_y)[1]) {
plot(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = la_mean_y[i,c(1,2,4,5:7,9,12)]/50, bg = alpha("green", 0.5), xlim = c(-350, 400), ylim = c(-350, 250), xaxt = "n", yaxt = "n")
points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ju_mean_y[i,c(1,2,4,5:7,9,12)]/25, bg = alpha("blue", 0.5))
points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ad_mean_y[i,c(1,2,4,5:7,9,12)]/5, bg = alpha("red", 0.5))
text(mds.coor[,1], -mds.coor[,2], rownames(mds.coor), font = 2)
abline(h=0,v=0, col = "grey", lty = 2)
legend("topright", rownames(la_mean_y)[i], bty = "n", text.font = 4)
}
#legend("bottomright", legend = c("Larvae (x10)", "Juveniles (x5)", "Adults"), fill = alpha(c("green", "blue", "red"), 0.5), bty = "n")
# And with a legend:
d <- as.dist(dist)
mds.coor <- cmdscale(d)
mds.coor <- mds.coor[c(1,2,4,5,6,7,9,13),]
#quartz("", 11, 7)
par(mar = c(0,0,0,0), mfrow = c(3,5))
for (i in 1:dim(ad_mean_y)[1]) {
plot(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = la_mean_y[i,c(1,2,4,5:7,9,12)]/50, bg = alpha("green", 0.5), xlim = c(-350, 400), ylim = c(-350, 210), xaxt = "n", yaxt = "n")
points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ju_mean_y[i,c(1,2,4,5:7,9,12)]/21, bg = alpha("blue", 0.5))
points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ad_mean_y[i,c(1,2,4,5:7,9,12)]/5, bg = alpha("red", 0.5))
a <- mds.coor[,1]+25
b <- -mds.coor[,2]+25
text(a,b, rownames(mds.coor), font = 2, cex=1.2)
abline(h=0,v=0, col = "grey", lty = 2)
legend("topright", rownames(la_mean_y)[i], bty = "n", text.font = 4)
legend("bottomright", legend = c("Tadpoles (x10)", "Subadults (x5)", "Adults", "Tadpoles (x10) & Adults", "Subadults (x5) & Adults", "Everything"), fill = alpha(c("green", "blue", "red", "saddlebrown", "mediumslateblue", "maroon4"), 0.5), bty = "n")
}
# Can I place the legend somewhere else?
d <- as.dist(dist)
mds.coor <- cmdscale(d)
mds.coor <- mds.coor[c(1,2,4,5,6,7,9,13),]
par(mar = c(0,0,0,0), mfrow = c(3,5))
for (i in 1:dim(ad_mean_y)[1]) {
plot(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = la_mean_y[i,c(1,2,4,5:7,9,12)]/50, bg = alpha("green", 0.5), xlim = c(-350, 400), ylim = c(-350, 210), xaxt = "n", yaxt = "n")
points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ju_mean_y[i,c(1,2,4,5:7,9,12)]/21, bg = alpha("blue", 0.5))
points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ad_mean_y[i,c(1,2,4,5:7,9,12)]/5, bg = alpha("red", 0.5))
a <- mds.coor[,1]+25
b <- -mds.coor[,2]+25
text(a,b, rownames(mds.coor), font = 2, cex=1.2)
abline(h=0,v=0, col = "grey", lty = 2)
legend("topright", inset=c(0,0), rownames(la_mean_y)[i], bty = "n", text.font = 4)
}
legend("bottomright", inset=c(0,0), legend = c("Tadpoles (x10)", "Subadults (x5)", "Adults", "T & A", "S & A", "All"), fill = alpha(c("green", "blue", "red", "saddlebrown", "mediumslateblue", "maroon4"), 0.5), bty = "n")
# A little bit more polishing:
d <- as.dist(dist)
mds.coor <- cmdscale(d)
mds.coor <- mds.coor[c(1,2,4,5,6,7,9,13),]
par(mar = c(0,0,0,0), mfrow = c(5,3))
for (i in 1:dim(ad_mean_y)[1]) {
plot(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = la_mean_y[i,c(1,2,4,5:7,9,12)]/50, bg = alpha("green", 0.5), xlim = c(-350, 450), ylim = c(-400, 270), xaxt = "n", yaxt = "n")
points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ju_mean_y[i,c(1,2,4,5:7,9,12)]/21, bg = alpha("blue", 0.5))
points(mds.coor[,1], -mds.coor[,2], type="p", xlab="", ylab="", pch = 21, cex = ad_mean_y[i,c(1,2,4,5:7,9,12)]/5, bg = alpha("red", 0.5))
a <- mds.coor[,1]+25
b <- -mds.coor[,2]+25
text(a,b, rownames(mds.coor), font = 2, cex=1.2)
abline(h=0,v=0, col = "grey", lty = 2)
legend("topright", inset=c(0,0), rownames(la_mean_y)[i], bty = "n", text.font = 4, cex = 1.5)
}
legend("bottomright", inset=c(0,0), legend = c("Tadpoles (x10)", "Subadults (x5)", "Adults", "T & A", "S & A", "All"), fill = alpha(c("green", "blue", "red", "saddlebrown", "mediumslateblue", "maroon4"), 0.5), bty = "n")
If I plot this last plot in R-Studio it looks quite OK. Legend only in the last plot and years clearly visible. In R Markdown it looks not so nice.
Now let’s make some useful summary plots. We decided to concentrate on the following life-stages: number of egg masses, tadpoles, subadults and adults.
I will plot tadpoles, subadults and adults in one plot and then egg masses in a separate plot. For tadpoles, subadults and adults I also show the log-transformed counts. For egg masses I don’t because there are so many zeros and the log of zero is an undefined number! –> Also: for eggs I only take counts during the breeding season!
Code for Supplementary Figure S4
# plotting with log
par(mfrow = c(1,3), oma=c(0,0,2,0))
boxplot(larvae ~ year, subset(ndat, larvae > 0), log = "y", ylab="", main = "Tadpoles (when present)", las = 2, cex.axis=0.75)
mtext("LOG(Mean counts per year)", side=2, line=3, cex=0.8)
boxplot(subadults ~ year, subset(ndat, subadults > 0), log = "y", main = "Subadults (when present)", las = 2)
boxplot(adults ~ year, subset(ndat, adults > 0), log = "y", main = "Adults (when present)", las = 2)
# trying to plot all years, without LOG, real counts, until 2015:
write.table(ndat, file="ndat_raw.txt")
# --> I manually added three years with zero counts in Excel!
ndat.v2 <- read.delim("~/DusyBasin/Input_summaryStats/ndat_2015.txt")
par(mfrow = c(1,3), oma=c(0,0,2,0))
boxplot(larvae ~ year, data=ndat.v2, las = 2, ylab="", main="Tadpoles", cex.axis=0.75)
mtext("Mean counts per year", side=2, line=3, cex=0.8)
boxplot(subadults ~ year, data=ndat.v2, main = "Subadults", las = 2)
boxplot(adults ~ year, data=ndat.v2, main = "Adults", las = 2)
Code for Supplementary Figures S5 and S7
I am plotting the eggs separately.
# for eggs take only breeding season:
data.breeding <- subset(ndat.v2, data$season=="b")
eg_mean_y_b <- tapply(data.breeding$eggs, list(data.breeding$year, data.breeding$site), mean, na.rm = TRUE)
ndat.egg <- data.frame(year = as.numeric(rep(rownames(eg_mean_y_b), times = dim(eg_mean_y_b)[2])), site = as.numeric(rep(colnames(eg_mean_y_b), each = dim(eg_mean_y_b)[1])), eggs = as.numeric(eg_mean_y_b))
summary(ndat.egg)## year site eggs
## Min. :2002 Min. : 1.000 Min. : 0.0000
## 1st Qu.:2003 1st Qu.: 4.000 1st Qu.: 0.0000
## Median :2006 Median : 7.500 Median : 0.0000
## Mean :2007 Mean : 9.214 Mean : 0.9900
## 3rd Qu.:2013 3rd Qu.:11.000 3rd Qu.: 0.1429
## Max. :2014 Max. :22.000 Max. :17.0000
## NA's :29# First a plot of the eggs only from the breeding season:
par(mfrow = c(1,2), oma=c(0,0,2,0))
boxplot(eggs ~ year, data=ndat.egg, las = 2, ylab="Counts", ylim=c(-1,15), main="Mean counts per year")
title("Mean number of egg masses per year", outer=TRUE)
boxplot(eggs ~ year, subset(ndat, eggs > 0), log = "y", ylab = "LOG(Counts)", las = 2, main="LOG(Mean counts per year when present)")
par(mfrow = c(1,2), oma=c(0,0,2,0))
boxplot(eggs ~ year, data=ndat, las = 2, ylab="Counts", ylim=c(-1,15), main="Mean counts per year")
title("Mean number of egg masses per year", outer=TRUE)
boxplot(eggs ~ year, subset(ndat, eggs > 0), log = "y", ylab = "LOG(Counts)", las = 2, main="LOG(Mean counts per year when present)")
Code for Supplementary Figure S6
Variation in counts among lakes at Dusy Basin. Again, I am first showing LOG-transformed values! Also: I am splitting up again and showing egg data only for breeding season!
par(mfrow = c(3,1), mar = c(3,3,2,1), oma = c(2,2,0,0))
boxplot(larvae ~ site, subset(ndat, larvae > 0), log = "y", main = "Tadpoles (when present)", las = 2)
boxplot(subadults ~ site, subset(ndat, subadults > 0), log = "y", main = "Subadults (when present)", las = 2)
boxplot(adults ~ site, subset(ndat, adults > 0), log = "y", main = "Adults (when present)", las = 1)
mtext("LOG(Mean counts per site)", 2, 0, TRUE)
mtext("Sites", 1, 0, TRUE)
# Now trying to plot this also for ALL lakes!
par(mfrow = c(1,3), mar = c(3,3,2,1), oma = c(2,2,0,0))
boxplot(larvae ~ site, data=ndat, main = "Tadpoles", las = 2)
boxplot(subadults ~ site, data=ndat, main = "Subadults", las = 2)
boxplot(adults ~ site, data=ndat, main = "Adults", las = 1)
mtext("Mean counts per site", 2, 0, TRUE)
mtext("Sites", 1, 0, TRUE)
# And now an individual plot for the eggs, again first only for the breeding season:
quartz("", 12, 6)
par(mfrow = c(1,2), mar = c(3,3,2,1), oma = c(2,2,0,0))
boxplot(eggs ~ site, data=ndat.egg, main = "Egg masses", las = 2)
boxplot(eggs ~ site, subset(ndat, eggs > 0), log = "y", main = "LOG(Egg masses)", las = 2)
mtext("Sites", 1, 0, TRUE, cex=1.2)
mtext("Counts", 2, 0, TRUE, cex=1.2)
# And then for all seasons:
quartz("", 12, 6)
par(mfrow = c(1,2), mar = c(3,3,2,1), oma = c(2,2,0,0))
boxplot(eggs ~ site, data=ndat, main = "Egg masses", las = 2)
boxplot(eggs ~ site, subset(ndat, eggs > 0), log = "y", main = "LOG(Egg masses when present)", las = 2)
mtext("Sites", 1, 0, TRUE, cex=1.2)
mtext("Counts", 2, 0, TRUE, cex=1.2)
Code for Supplementary Figure S8
How does fish presence affect frog abundance at different life stages?
# Lake breeding quality
eg_max_y <- tapply(data$eggs, list(data$year, data$site), max, na.rm = TRUE)
breed_eggs <- apply(eg_max_y, 2, mean, na.rm = TRUE) / sum(apply(eg_max_y, 2, mean, na.rm = TRUE))
la_max_y <- tapply(data$larvae, list(data$year, data$site), max, na.rm = TRUE)
breed_larv <- apply(la_max_y, 2, mean, na.rm = TRUE) / sum(apply(la_max_y, 2, mean, na.rm = TRUE))
#quartz("", 12, 5)
par(mfrow = c(1,2), mar = c(5,4,1,1))
barplot(rbind(breed_eggs, breed_larv), beside = TRUE, las = 2)
mtext("Average number of egg masses / larvae (standardized)", 2, 3, las = 0)
plot(breed_larv[-c(9, 11, 13)] ~ apply(min_surf, 2, mean, na.rm = TRUE), pch = 16, log = "x", xlab = "Lake surface", xlim = c(0.3, 150000))
points(breed_eggs[-c(9,11,13)] ~ apply(min_surf, 2, mean, na.rm = TRUE), pch = 17, log = "x")## Warning in plot.xy(xy.coords(x, y), type = type, ...): "log" is not a
## graphical parametertext(apply(min_surf, 2, mean, na.rm = TRUE), breed_larv[-c(9,11,13)], paste("Lake n°", names(breed_eggs[-c(9,11,13)])), pos = c(3, 4, 4, 4, 4, 4,3 , 4,3 ,4 , 4))
# Hypothesis 1: Best lakes are larger lakes / but fish eat frogs
fish <- table(data$fish, data$site, data$year)
test <- fish[2,,] / (fish[1,,] + fish[2,,])
test <- apply(test, 1, mean, na.rm = TRUE)
fish2 <- table(data$fish, data$site)
test2 <- fish2[2,] / (fish2[1,] + fish2[2,])
par(mfrow = c(1,2), mar = c(5,4,1,1), las = 1)
boxplot(min_surf+1, log = "y", las = 2, ylab = "Minimal lake surface")
#barplot(t(test), beside = TRUE, ylab = "Proportion of observed fish")
barplot(rbind(test, test2), beside = TRUE, ylab = "Proportion of observed fish", las =2)
lake_quality <- data.frame(breed_larv, breed_eggs, lake = as.numeric(names(breed_larv)))
prop_fish <- data.frame(prop_fish = test, lake = as.numeric(names(test)))
surface <- data.frame(surf.min = apply(min_surf, 2, mean, na.rm = TRUE), surf.mean = apply(mean_surf, 2, mean, na.rm = TRUE), surf.max = apply(max_surf, 2, max, na.rm = TRUE), lake = as.numeric(colnames(mean_surf)))
lake_quality <- merge(lake_quality, prop_fish)
lake_quality <- merge(lake_quality, surface)
lake_quality <- lake_quality[order(lake_quality$lake),]
lake_quality$fish <- "fishless"
lake_quality$fish[lake_quality$lake %in% c(1,3)] <- "fishes"
# Lake 1, lake 3, early (lake 20 - fish died in 2000)
par(las = 1, plt = c(0,0.7, 0, 1), oma = c(4,5,1,1))
plot(breed_larv ~ surf.max, subset(lake_quality, fish == "fishless"), pch = 16, ylab = "", xlab = "", ylim = c(0, 0.5), xlim = c(0, 2800))
points(breed_eggs ~ surf.max, lake_quality, pch = 17)
res <- lm(breed_larv ~ surf.max + I(surf.max^2), subset(lake_quality, fish == "fishless"))
res <- lm(breed_larv ~ I(surf.max^2), subset(lake_quality, fish == "fishless"))
pred <- predict(res, newdata = data.frame(surf.max = c(0:3000)))
points(pred ~ c(0:3000), type = "l")
text(lake_quality$surf.max, lake_quality$breed_larv, lake_quality$lake, pos = c(rep(3, 7), 4, rep(3, 3)))
points(breed_larv ~ surf.max, subset(lake_quality, fish == "fishes"), pch = 1)
par(las = 1, plt = c(0.72,1, 0, 1), new = TRUE)
plot(breed_larv ~ surf.max, lake_quality, pch = 1, ylab = "", xlab = "", ylim = c(0, 0.5), xlim = c(50500, 51500), yaxt = "n", xaxt = "n")
axis(1, at = c(50500, 51000, 51500))
text(lake_quality$surf.max, lake_quality$breed_larv, lake_quality$lake, pos = 3)
mtext("Average number of tadpoles (standardized)", 2, 3, TRUE, las = 0)
mtext(expression(paste("Lake surface [", m^2, "]", sep = "")), 1, 3, TRUE)
# And as Steve Beissinger suggested let's try to also plot number of tadpoles in logarithm:
par(las = 1, plt = c(0,0.7, 0, 1), oma = c(4,5,1,1))
breed_larv <- log(breed_larv)
plot(breed_larv ~ surf.max, subset(lake_quality, fish == "fishless"), pch = 16, ylab = "", xlab = "", ylim = c(0, 0.5), xlim = c(0, 2800))
# points(breed_eggs ~ surf.max, lake_quality, pch = 17)
res <- lm(breed_larv ~ surf.max + I(surf.max^2), subset(lake_quality, fish == "fishless"))
res <- lm(breed_larv ~ I(surf.max^2), subset(lake_quality, fish == "fishless"))
pred <- predict(res, newdata = data.frame(surf.max = c(0:3000)))
points(pred ~ c(0:3000), type = "l")
text(lake_quality$surf.max, lake_quality$breed_larv, lake_quality$lake, pos = c(rep(3, 7), 4, rep(3, 3)))
points(breed_larv ~ surf.max, subset(lake_quality, fish == "fishes"), pch = 1)
par(las = 1, plt = c(0.72,1, 0, 1), new = TRUE)
plot(breed_larv ~ surf.max, lake_quality, pch = 1, ylab = "", xlab = "", ylim = c(0, 0.5), xlim = c(50500, 51500), yaxt = "n", xaxt = "n")
axis(1, at = c(50500, 51000, 51500))
text(lake_quality$surf.max, lake_quality$breed_larv, lake_quality$lake, pos = 3)
mtext("Average number of tadpoles (standardized)", 2, 3, TRUE, las = 0)
mtext(expression(paste("Lake surface [", m^2, "]", sep = "")), 1, 3, TRUE)
From here on, I am not running this code for publication anymore. It was used for data exploration. However, if you want to see all the analyses we did, turn eval=TRUE and it will run like a charm.
Now we are moving on to environmental factors!
Hypothesis: The years with low snow cover, or drought, are when lakes have dried and when tadpoles have died.
Investigate for all different life stages:
wc <- tapply(dsnow$watcont, dsnow$year3, max, na.rm = TRUE)
wc <- data.frame(year = as.numeric(rownames(wc)), wc.max = wc, wc.mean = tapply(dsnow$watcont, dsnow$year3, mean, na.rm = TRUE))
# Season length estimated from figure with daily data
# ndat = new data, max per site and year
no_data2004 <- ndat[1:14,]
no_data2004$site <- unique(ndat$site)
no_data2004$year <- 2004
no_data2004[,3:9] <- NA
no_data1996 <- no_data2004; no_data1996$year <- 1996
no_data1995 <- no_data2004; no_data1995$year <- 1995
no_data1994 <- no_data2004; no_data1994$year <- 1994
no_data1993 <- no_data2004; no_data1993$year <- 1993
no_data1992 <- no_data2004; no_data1992$year <- 1992
ndat <- rbind(ndat, no_data1992, no_data1993, no_data1994, no_data1995, no_data1996, no_data2004)
data2 <- merge(ndat, sle, all = TRUE)
data2 <- merge(data2, wc, all = TRUE)
# What environmental factor explains eggs abundance
res <- step(lmer(eggs ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 13 -341.78 709.57
## (1 | site) 0 12 -359.35 742.71 35.145 1 3.061e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## sum.years:wc.max 1 0.133 0.133 1 128.42 0.0208
## sum.years:win.years 2 2.253 2.253 1 129.41 0.3567
## win.years:wc.max 3 0.163 0.163 1 130.69 0.0260
## win.years:winter 4 18.239 18.239 1 131.01 2.9243
## sum.years:winter 5 6.844 6.844 1 132.34 1.0824
## winter:wc.max 6 21.473 21.473 1 133.53 3.3928
## wc.max 7 12.900 12.900 1 134.91 1.9992
## winter 8 9.120 9.120 1 135.07 1.4002
## sum.years 0 115.062 115.062 1 136.02 17.6067
## win.years 0 53.193 53.193 1 136.04 8.1396
## Pr(>F)
## sum.years:wc.max 0.885482
## sum.years:win.years 0.551392
## win.years:wc.max 0.872230
## win.years:winter 0.089622 .
## sum.years:winter 0.300056
## winter:wc.max 0.067702 .
## wc.max 0.159687
## winter 0.238764
## sum.years 4.873e-05 ***
## win.years 0.005008 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## eggs ~ sum.years + win.years + (1 | site)res <- lmer(eggs ~ sum.years + win.years + ( 1|site), subset(data2, year < 2010)); summary(res); anova(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: eggs ~ sum.years + win.years + (1 | site)
## Data: subset(data2, year < 2010)
##
## REML criterion at convergence: 731
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3686 -0.2727 -0.0549 0.1046 6.8095
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 3.340 1.828
## Residual 6.535 2.556
## Number of obs: 152, groups: site, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.465 3.531 141.600 4.096 7.03e-05 ***
## sum.years -15.630 3.725 136.021 -4.196 4.87e-05 ***
## win.years -11.092 3.888 136.036 -2.853 0.00501 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sm.yrs
## sum.years -0.904
## win.years -0.911 0.685## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## sum.years 115.062 115.062 1 136.02 17.6067 4.873e-05 ***
## win.years 53.193 53.193 1 136.04 8.1396 0.005008 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1a <- 100
pred <- predict(res, newdata = data.frame(sum.years = rep(seq(0.2, 0.8, length.out = a), each = a), win.years = rep(seq(0.2, 0.8, length.out = a), times = a), site = rep(2, a*a)))
par(mar = c(4,4,1,1), las = 1)
image(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), xlab = "Summer length", ylab = "Winter length")
contour(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), add = TRUE, drawlabels = TRUE)
# What environmental factor explains larvae abundance
res <- step(lmer(larvae ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 13 -963.25 1952.5
## (1 | site) 0 12 -992.64 2009.3 58.768 1 1.774e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## sum.years:winter 1 12 12 1 128.38 0.0003 0.98624
## win.years:winter 2 1346 1346 1 129.10 0.0329 0.85645
## winter:wc.max 3 9065 9065 1 130.37 0.2229 0.63761
## win.years:wc.max 4 13515 13515 1 131.36 0.3343 0.56410
## sum.years:win.years 5 5733 5733 1 132.71 0.1425 0.70640
## winter 6 144879 144879 1 133.13 3.6244 0.05910
## win.years 0 224560 224560 1 134.21 5.5112 0.02036
## sum.years:wc.max 0 215705 215705 1 134.32 5.2939 0.02294
##
## sum.years:winter
## win.years:winter
## winter:wc.max
## win.years:wc.max
## sum.years:win.years
## winter .
## win.years *
## sum.years:wc.max *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## larvae ~ sum.years + win.years + wc.max + (1 | site) + sum.years:wc.maxres <- lmer(larvae ~ sum.years + win.years + ( 1|site), subset(data2, year < 2010)); summary(res); anova(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: larvae ~ sum.years + win.years + (1 | site)
## Data: subset(data2, year < 2010)
##
## REML criterion at convergence: 2042.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7947 -0.2347 -0.1101 0.1222 7.3332
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 34916 186.9
## Residual 41878 204.6
## Number of obs: 152, groups: site, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 670.7 284.4 144.0 2.359 0.0197 *
## sum.years -445.2 298.2 136.1 -1.493 0.1377
## win.years -638.5 311.3 136.1 -2.051 0.0422 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sm.yrs
## sum.years -0.898
## win.years -0.906 0.685## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## sum.years 93354 93354 1 136.13 2.2292 0.13774
## win.years 176226 176226 1 136.14 4.2080 0.04215 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pred <- predict(res, newdata = data.frame(sum.years = rep(seq(0.2, 0.8, length.out = a), each = a), win.years = rep(seq(0.2, 0.8, length.out = a), times = a), site = rep(2, a*a)))
#quartz()
par(mar = c(4,4,1,1), las = 1)
image(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), xlab = "Summer length", ylab = "Winter length")
contour(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), add = TRUE, drawlabels = TRUE)
res <- lmer(larvae ~ sum.years*wc.max + ( 1|site), subset(data2, year < 2010)); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: larvae ~ sum.years * wc.max + (1 | site)
## Data: subset(data2, year < 2010)
##
## REML criterion at convergence: 2044.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8225 -0.2484 -0.1318 0.1175 7.1305
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 36429 190.9
## Residual 42066 205.1
## Number of obs: 152, groups: site, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 752.40 500.29 138.29 1.504 0.1349
## sum.years -1355.69 919.61 135.36 -1.474 0.1427
## wc.max -18.30 11.56 135.37 -1.583 0.1158
## sum.years:wc.max 41.71 22.80 135.34 1.829 0.0696 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sm.yrs wc.max
## sum.years -0.983
## wc.max -0.945 0.969
## sm.yrs:wc.m 0.869 -0.927 -0.973pred <- predict(res, newdata = data.frame(sum.years = rep(seq(0.2, 0.8, length.out = a), each = a), wc.max = rep(seq(10, 60, length.out = a), times = a), site = rep(2, a*a)))
image(matrix(pred, nrow = a), xlab = "Summer length", ylab = "Water max")
contour(matrix(pred, nrow = a), add = TRUE, drawlabels = TRUE)
# What environmental factor explains juvenile abundance
res <- step(lmer(juveniles ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 13 -756.14 1538.3
## (1 | site) 0 12 -777.45 1578.9 42.611 1 6.68e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## win.years:winter 1 53 53 1 128.26 0.0239
## sum.years:wc.max 2 1496 1496 1 128.96 0.6744
## winter:wc.max 3 2557 2557 1 130.11 1.1547
## sum.years:win.years 0 44577 44577 1 131.89 20.1195
## sum.years:winter 0 18088 18088 1 131.35 8.1639
## win.years:wc.max 0 31898 31898 1 131.42 14.3967
## Pr(>F)
## win.years:winter 0.877434
## sum.years:wc.max 0.413051
## winter:wc.max 0.284552
## sum.years:win.years 1.566e-05 ***
## sum.years:winter 0.004970 **
## win.years:wc.max 0.000225 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## juveniles ~ sum.years + win.years + winter + wc.max + (1 | site) +
## sum.years:win.years + sum.years:winter + win.years:wc.max res <- lmer(juveniles ~ sum.years * win.years + ( 1|site), subset(data2, year < 2010)); summary(res); anova(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: juveniles ~ sum.years * win.years + (1 | site)
## Data: subset(data2, year < 2010)
##
## REML criterion at convergence: 1604
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6173 -0.2714 -0.1482 0.1142 5.8718
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 1575 39.69
## Residual 2513 50.13
## Number of obs: 152, groups: site, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 320.0 234.0 136.5 1.368 0.174
## sum.years -470.4 424.9 136.0 -1.107 0.270
## win.years -563.6 419.8 136.0 -1.342 0.182
## sum.years:win.years 897.7 780.5 136.0 1.150 0.252
##
## Correlation of Fixed Effects:
## (Intr) sm.yrs wn.yrs
## sum.years -0.987
## win.years -0.988 0.990
## sm.yrs:wn.y 0.955 -0.985 -0.983## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## sum.years 3080.3 3080.3 1 135.98 1.2255 0.2702
## win.years 4528.8 4528.8 1 136.03 1.8018 0.1817
## sum.years:win.years 3325.4 3325.4 1 136.05 1.3230 0.2521pred <- predict(res, newdata = data.frame(sum.years = rep(seq(0.2, 0.8, length.out = a), each = a), win.years = rep(seq(0.2, 0.8, length.out = a), times = a), site = rep(2, a*a)))
par(mar = c(4,4,1,1), las = 1)
image(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), xlab = "Summer length", ylab = "Winter length")
contour(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), add = TRUE, drawlabels = TRUE)
# What environmental factor explains subadults abundance
res <- step(lmer(subadults ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 13 -578.76 1183.5
## (1 | site) 0 12 -588.81 1201.6 20.104 1 7.336e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## sum.years:win.years 0 1418.48 1418.48 1 129.00 7.4538
## sum.years:winter 0 896.33 896.33 1 128.65 4.7100
## sum.years:wc.max 0 788.27 788.27 1 128.60 4.1422
## win.years:winter 0 874.24 874.24 1 128.64 4.5939
## win.years:wc.max 0 1256.32 1256.32 1 128.76 6.6017
## winter:wc.max 0 1704.74 1704.74 1 128.12 8.9580
## Pr(>F)
## sum.years:win.years 0.007216 **
## sum.years:winter 0.031824 *
## sum.years:wc.max 0.043882 *
## win.years:winter 0.033967 *
## win.years:wc.max 0.011328 *
## winter:wc.max 0.003315 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## subadults ~ (sum.years + win.years + winter + wc.max)^2 + (1 |
## site) res <- lmer(subadults ~ sum.years * win.years + ( 1|site), subset(data2, year < 2010)); summary(res); anova(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: subadults ~ sum.years * win.years + (1 | site)
## Data: subset(data2, year < 2010)
##
## REML criterion at convergence: 1231.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4208 -0.5041 -0.0731 0.0982 6.3068
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 69.41 8.331
## Residual 212.85 14.589
## Number of obs: 152, groups: site, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -102.88 67.98 136.87 -1.513 0.1325
## sum.years 219.77 123.47 136.63 1.780 0.0773 .
## win.years 180.60 122.00 136.71 1.480 0.1411
## sum.years:win.years -362.92 226.80 136.73 -1.600 0.1119
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sm.yrs wn.yrs
## sum.years -0.987
## win.years -0.989 0.990
## sm.yrs:wn.y 0.955 -0.985 -0.983## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## sum.years 674.32 674.32 1 136.63 3.1680 0.07732 .
## win.years 466.44 466.44 1 136.71 2.1914 0.14109
## sum.years:win.years 545.04 545.04 1 136.73 2.5607 0.11186
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pred <- predict(res, newdata = data.frame(sum.years = rep(seq(0.2, 0.8, length.out = a), each = a), win.years = rep(seq(0.2, 0.8, length.out = a), times = a), site = rep(2, a*a)))
par(mar = c(4,4,1,1), las = 1)
image(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), xlab = "Summer length", ylab = "Winter length")
contour(seq(0.2, 0.8, length.out = a), seq(0.2, 0.8, length.out = a), matrix(pred, nrow = a), add = TRUE, drawlabels = TRUE)
# What environmental factor explains adults abundance
res <- step(lmer(adults ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 13 -544.36 1114.7
## (1 | site) 0 12 -588.25 1200.5 87.781 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## winter:wc.max 1 102.48 102.48 1 128.02 0.9872
## sum.years:winter 2 36.16 36.16 1 129.21 0.3485
## sum.years:wc.max 3 14.40 14.40 1 130.12 0.1395
## sum.years:win.years 4 68.68 68.68 1 131.33 0.6698
## sum.years 5 15.99 15.99 1 132.24 0.1564
## win.years:winter 0 473.20 473.20 1 133.02 4.6599
## win.years:wc.max 0 1830.34 1830.34 1 133.23 18.0244
## Pr(>F)
## winter:wc.max 0.32230
## sum.years:winter 0.55601
## sum.years:wc.max 0.70936
## sum.years:win.years 0.41459
## sum.years 0.69315
## win.years:winter 0.03267 *
## win.years:wc.max 4.062e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## adults ~ win.years + winter + wc.max + (1 | site) + win.years:winter +
## win.years:wc.max# What environmental factor explains adults and subadults abundance
res <- step(lmer(frogs ~ (sum.years + win.years + winter + wc.max)^2 + (1|site), subset(data2, year < 2010)));res## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 13 -640.94 1307.9
## (1 | site) 0 12 -671.15 1366.3 60.424 1 7.648e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## sum.years:wc.max 1 1555.6 1555.6 1 128.20 3.6733
## sum.years:winter 2 325.6 325.6 1 128.97 0.7545
## win.years:winter 3 270.1 270.1 1 129.97 0.6274
## sum.years:win.years 0 1890.9 1890.9 1 131.34 4.4069
## win.years:wc.max 0 5852.7 5852.7 1 131.09 13.6404
## winter:wc.max 0 1873.0 1873.0 1 131.37 4.3653
## Pr(>F)
## sum.years:wc.max 0.0575177 .
## sum.years:winter 0.3866705
## win.years:winter 0.4297400
## sum.years:win.years 0.0377068 *
## win.years:wc.max 0.0003238 ***
## winter:wc.max 0.0386071 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## frogs ~ sum.years + win.years + winter + wc.max + (1 | site) +
## sum.years:win.years + win.years:wc.max + winter:wc.maxDrought and season length investigations
How season length affects drying out of pools and lakes:
table(dry$dryness, dry$lake)##
## 2 5 6 7 8 9 10 11
## dry 8 0 8 0 11 7 10 7
## low 2 3 0 5 0 1 0 0
## wet 4 11 6 9 3 6 4 2test <- merge(sle, dry)
res <- lmer(winter ~ dryness + (1|lake), test); anova(res)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## dryness 0.070291 0.035146 2 104 4.8678 0.009531 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res <- lmer(sum.years ~ dryness + (1|lake), test); anova(res)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## dryness 0.05313 0.026565 2 104 4.2857 0.01627 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1par(mfrow = c(1,2), las = 1, mar = c(3,4,1,1), oma = c(0,0,2,1))
boxplot(winter ~ dryness, subset(test, lake == 2), ylab = "Winter length")
res <- lm(winter ~ dryness, subset(test, lake == 2)); anova(res)## Analysis of Variance Table
##
## Response: winter
## Df Sum Sq Mean Sq F value Pr(>F)
## dryness 2 0.050474 0.0252369 4.9102 0.02992 *
## Residuals 11 0.056537 0.0051397
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 res2 <- lmer(winter ~ dryness + (1|lake), test); anova(res2)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## dryness 0.070291 0.035146 2 104 4.8678 0.009531 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 legend("topleft", c(paste(p.trans(anova(res)$Pr[1]), "(lake 2)"), paste(p.trans(anova(res2)$Pr[1]), "(all lakes)")), bty = "n")
boxplot(sum.years ~ dryness, subset(test, lake == 2), ylab = "Summer length")
res <- lm(sum.years ~ dryness, subset(test, lake == 2)); anova(res)## Analysis of Variance Table
##
## Response: sum.years
## Df Sum Sq Mean Sq F value Pr(>F)
## dryness 2 0.024695 0.0123477 2.0327 0.1773
## Residuals 11 0.066821 0.0060746 res2 <- lmer(sum.years ~ dryness + (1|lake), test); anova(res2)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## dryness 0.05313 0.026565 2 104 4.2857 0.01627 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 legend("bottomleft", c(paste(p.trans(anova(res)$Pr[1]), "(lake 2)"), paste(p.trans(anova(res2)$Pr[1]), "(all lakes)")), bty = "n")
mtext("Lake 2 example", 3, 0, font = 2, out = TRUE)
Relationship between drought (in drying lakes) and snow pack and season lengths.
# Combined files year of drought versus year after drought
test <- merge(sle, dry)
drought <- merge(dry, data)
drought <- subset(drought, season == "b")
dry2 <- test
dry2$year <- I(dry2$year + 1)
drought2 <- merge(dry2, data)
drought2 <- subset(drought2, season == "b")
dry$site <- dry$lake
drought <- merge(dry, data2)
dry2 <- dry
dry2$year <- I(dry$year + 1)
drought2 <- merge(dry2, data2)
# Data structure
table(drought$year, drought$site)##
## 2 5 6 7 8 9 10 11
## 1997 1 1 1 1 1 1 1 1
## 1998 1 1 1 1 1 1 1 1
## 1999 1 1 1 1 1 1 1 1
## 2000 1 1 1 1 1 1 1 1
## 2001 1 1 1 1 1 1 1 1
## 2002 1 1 1 1 1 1 1 1
## 2003 1 1 1 1 1 1 1 1
## 2004 1 1 1 1 1 1 1 1
## 2005 1 1 1 1 1 1 1 1
## 2006 1 1 1 1 1 1 1 1
## 2007 1 1 1 1 1 1 1 1
## 2008 1 1 1 1 1 1 1 1
## 2009 1 1 1 1 1 1 1 1
## 2010 1 1 1 1 1 1 1 1table(drought$site, drought$dryness)##
## dry low wet
## 2 8 2 4
## 5 0 3 11
## 6 8 0 6
## 7 0 5 9
## 8 11 0 3
## 9 7 1 6
## 10 10 0 4
## 11 7 0 2# Combined files year of drought versus year after drought
par(mfrow = c(2,4), mar = c(1,2,2,1), oma = c(4,3,1,0))
boxplot(I(eggs+1) ~ dryness, drought, log = "y", main = "Egg masses", las = 1)
res <- lmer(log(eggs+1) ~ dryness + (1|site), drought); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(eggs + 1) ~ dryness + (1 | site)
## Data: drought
##
## REML criterion at convergence: 175.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.95636 -0.27087 -0.06274 -0.02951 2.80811
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.1438 0.3793
## Residual 0.3175 0.5634
## Number of obs: 94, groups: site, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.28249 0.16258 9.76009 1.738 0.114
## drynesslow 0.07577 0.25000 90.98213 0.303 0.763
## drynesswet 0.11480 0.14308 90.99677 0.802 0.424
##
## Correlation of Fixed Effects:
## (Intr) drynssl
## drynesslow -0.284
## drynesswet -0.421 0.433 legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Pop. size (during drought)", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ dryness, drought, log = "y", main = "Larvae", las = 1)
res <- lmer(log(larvae +1) ~ dryness + (1|site), drought); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(larvae + 1) ~ dryness + (1 | site)
## Data: drought
##
## REML criterion at convergence: 382
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.03642 -0.57303 -0.05086 0.57345 2.43186
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 3.666 1.915
## Residual 2.851 1.688
## Number of obs: 94, groups: site, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.1603 0.7325 8.0076 2.949 0.0184 *
## drynesslow 0.7316 0.7663 88.2544 0.955 0.3423
## drynesswet 0.1942 0.4393 88.4867 0.442 0.6595
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) drynssl
## drynesslow -0.198
## drynesswet -0.288 0.456 legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(juveniles+1) ~ dryness, drought, log = "y", main = "Juveniles", las = 1)
res <- lmer(log(juveniles +1) ~ dryness + (1|site), drought); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(juveniles + 1) ~ dryness + (1 | site)
## Data: drought
##
## REML criterion at convergence: 305
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6818 -0.3661 -0.1764 0.1905 2.7698
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 1.016 1.008
## Residual 1.264 1.124
## Number of obs: 94, groups: site, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.72073 0.40165 8.50316 1.794 0.108
## drynesslow 0.27373 0.50614 89.75381 0.541 0.590
## drynesswet 0.04649 0.28998 89.99626 0.160 0.873
##
## Correlation of Fixed Effects:
## (Intr) drynssl
## drynesslow -0.236
## drynesswet -0.346 0.448 legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(frogs+1) ~ dryness, drought, log = "y", main = "Adults and subadults", las = 1)
res <- lmer(log(frogs +1) ~ dryness + (1|site), drought); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(frogs + 1) ~ dryness + (1 | site)
## Data: drought
##
## REML criterion at convergence: 267.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.10960 -0.51168 -0.02874 0.57383 2.39248
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.5294 0.7276
## Residual 0.8540 0.9241
## Number of obs: 94, groups: site, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.3780 0.2986 9.0854 4.614 0.00123 **
## drynesslow 0.4244 0.4136 90.4971 1.026 0.30760
## drynesswet 0.4616 0.2368 90.6815 1.949 0.05438 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) drynssl
## drynesslow -0.258
## drynesswet -0.380 0.442 legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(eggs+1) ~ dryness, drought2, log = "y", main = "", las = 1)
res <- lmer(log(eggs+1) ~ dryness + (1|site), drought2); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(eggs + 1) ~ dryness + (1 | site)
## Data: drought2
##
## REML criterion at convergence: 172.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.11230 -0.42444 0.01234 0.11625 2.83040
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.1360 0.3688
## Residual 0.3115 0.5581
## Number of obs: 93, groups: site, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.2013 0.1602 10.5685 1.256 0.2361
## drynesslow 0.2518 0.2195 89.9930 1.147 0.2543
## drynesswet 0.2625 0.1399 89.7020 1.876 0.0639 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) drynssl
## drynesslow -0.317
## drynesswet -0.430 0.453 legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Pop. size (after drought)", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ dryness, drought2, log = "y", main = "", las = 1)
res <- lmer(log(larvae +1) ~ dryness + (1|site), drought2); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(larvae + 1) ~ dryness + (1 | site)
## Data: drought2
##
## REML criterion at convergence: 386.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.95250 -0.59904 -0.04415 0.63750 2.32732
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 3.246 1.802
## Residual 3.166 1.779
## Number of obs: 93, groups: site, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.0328 0.7043 8.5488 2.886 0.019 *
## drynesslow 1.1589 0.7140 87.9479 1.623 0.108
## drynesswet 0.2407 0.4535 87.2611 0.531 0.597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) drynssl
## drynesslow -0.238
## drynesswet -0.318 0.471 legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(juveniles+1) ~ dryness, drought2, log = "y", main = "", las = 1)
res <- lmer(log(juveniles +1) ~ dryness + (1|site), drought2); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(juveniles + 1) ~ dryness + (1 | site)
## Data: drought2
##
## REML criterion at convergence: 300.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6295 -0.3508 -0.1245 0.1737 2.8947
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.881 0.9386
## Residual 1.253 1.1196
## Number of obs: 93, groups: site, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.5977 0.3815 9.0381 1.567 0.151
## drynesslow 0.6498 0.4459 89.1122 1.457 0.149
## drynesswet 0.2507 0.2836 88.4128 0.884 0.379
##
## Correlation of Fixed Effects:
## (Intr) drynssl
## drynesslow -0.273
## drynesswet -0.367 0.464 legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(frogs+1) ~ dryness, drought2, log = "y", main = "", las = 1)
res <- lmer(log(frogs +1) ~ dryness + (1|site), drought2); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(frogs + 1) ~ dryness + (1 | site)
## Data: drought2
##
## REML criterion at convergence: 272.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.88822 -0.43647 -0.03915 0.60032 2.47827
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.4261 0.6528
## Residual 0.9492 0.9743
## Number of obs: 93, groups: site, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.1324 0.2824 10.2815 4.011 0.00234 **
## drynesslow 0.7393 0.3835 89.9802 1.928 0.05705 .
## drynesswet 0.5953 0.2444 89.6422 2.436 0.01684 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) drynssl
## drynesslow -0.315
## drynesswet -0.426 0.454 legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Dry versus not dry lake", side = 1, line = 2, outer = TRUE)
Same figure for lake 2 only
par(mfrow = c(2,4), mar = c(1,2,2,1), oma = c(4,3,1,0))
boxplot(I(eggs+1) ~ dryness, subset(drought, site == 2), log = "y", main = "Egg masses", las = 1)
res <- lm(log(eggs+1) ~ dryness, subset(drought, site == 2)); summary(res)##
## Call:
## lm(formula = log(eggs + 1) ~ dryness, data = subset(drought,
## site == 2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2184 -1.2043 -0.2740 0.9563 1.3465
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.21844 0.44934 2.712 0.0219 *
## drynesslow -0.38932 0.95320 -0.408 0.6916
## drynesswet -0.01419 0.74515 -0.019 0.9852
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.189 on 10 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.0174, Adjusted R-squared: -0.1791
## F-statistic: 0.08854 on 2 and 10 DF, p-value: 0.916 legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
mtext("Pop. size (during subset(drought, site == 2))", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ dryness, subset(drought, site == 2), log = "y", main = "Larvae", las = 1)
res <- lm(log(larvae +1) ~ dryness, subset(drought, site == 2)); summary(res)##
## Call:
## lm(formula = log(larvae + 1) ~ dryness, data = subset(drought,
## site == 2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.635 -2.820 1.220 1.964 3.016
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.635 1.146 4.044 0.00235 **
## drynesslow -1.619 2.431 -0.666 0.52056
## drynesswet 1.026 1.901 0.540 0.60127
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.032 on 10 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.09267, Adjusted R-squared: -0.08879
## F-statistic: 0.5107 on 2 and 10 DF, p-value: 0.6149 legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
boxplot(I(juveniles+1) ~ dryness, subset(drought, site == 2), log = "y", main = "Juveniles", las = 1)
res <- lm(log(juveniles +1) ~ dryness, subset(drought, site == 2)); summary(res)##
## Call:
## lm(formula = log(juveniles + 1) ~ dryness, data = subset(drought,
## site == 2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0636 -0.6287 -0.3102 0.8974 2.0636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6287 0.4526 1.389 0.195
## drynesslow 1.4349 0.9602 1.494 0.166
## drynesswet -0.1685 0.7506 -0.225 0.827
##
## Residual standard error: 1.198 on 10 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2121, Adjusted R-squared: 0.05453
## F-statistic: 1.346 on 2 and 10 DF, p-value: 0.3036 legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
boxplot(I(frogs+1) ~ dryness, subset(drought, site == 2), log = "y", main = "Adults and subadults", las = 1)
res <- lm(log(frogs +1) ~ dryness, subset(drought, site == 2)); summary(res)##
## Call:
## lm(formula = log(frogs + 1) ~ dryness, data = subset(drought,
## site == 2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5423 -0.6908 0.4268 0.8693 2.5423
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.48253 0.58094 4.273 0.00163 **
## drynesslow 0.05973 1.23237 0.048 0.96230
## drynesswet 0.96780 0.96339 1.005 0.33879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.537 on 10 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.09669, Adjusted R-squared: -0.08398
## F-statistic: 0.5352 on 2 and 10 DF, p-value: 0.6014 legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
boxplot(I(eggs+1) ~ dryness, subset(drought2, site == 2), log = "y", main = "", las = 1)
res <- lm(log(eggs+1) ~ dryness, subset(drought2, site == 2)); summary(res)##
## Call:
## lm(formula = log(eggs + 1) ~ dryness, data = subset(drought2,
## site == 2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.28247 -0.76013 -0.06698 0.42889 1.50855
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7601 0.3591 2.117 0.0603 .
## drynesslow 0.5223 0.8029 0.651 0.5300
## drynesswet 1.3593 0.6875 1.977 0.0762 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.016 on 10 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.283, Adjusted R-squared: 0.1396
## F-statistic: 1.973 on 2 and 10 DF, p-value: 0.1895 legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
mtext("Pop. size (after subset(drought, site == 2))", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ dryness, subset(drought2, site == 2), log = "y", main = "", las = 1)
res <- lm(log(larvae +1) ~ dryness, subset(drought2, site == 2)); summary(res)##
## Call:
## lm(formula = log(larvae + 1) ~ dryness, data = subset(drought2,
## site == 2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0253 -3.2673 0.8187 1.9769 3.5987
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.0253 1.1312 3.558 0.0052 **
## drynesslow -0.7581 2.5295 -0.300 0.7705
## drynesswet 2.4870 2.1661 1.148 0.2776
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.2 on 10 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1423, Adjusted R-squared: -0.02922
## F-statistic: 0.8296 on 2 and 10 DF, p-value: 0.4641 legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
boxplot(I(juveniles+1) ~ dryness, subset(drought2, site == 2), log = "y", main = "", las = 1)
res <- lm(log(juveniles +1) ~ dryness, subset(drought2, site == 2)); summary(res)##
## Call:
## lm(formula = log(juveniles + 1) ~ dryness, data = subset(drought2,
## site == 2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.07709 -0.51713 -0.06515 0.02704 2.05004
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5171 0.3773 1.371 0.2005
## drynesslow -0.5171 0.8437 -0.613 0.5536
## drynesswet 1.5600 0.7225 2.159 0.0562 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.067 on 10 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.3743, Adjusted R-squared: 0.2492
## F-statistic: 2.991 on 2 and 10 DF, p-value: 0.09589 legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
boxplot(I(frogs+1) ~ dryness, subset(drought2, site == 2), log = "y", main = "", las = 1)
res <- lm(log(frogs +1) ~ dryness, subset(drought2, site == 2)); summary(res)##
## Call:
## lm(formula = log(frogs + 1) ~ dryness, data = subset(drought2,
## site == 2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1841 -0.6057 0.4550 1.1636 1.1763
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1841 0.4826 4.526 0.0011 **
## drynesslow -1.0205 1.0791 -0.946 0.3666
## drynesswet 1.9939 0.9241 2.158 0.0563 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.365 on 10 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.4075, Adjusted R-squared: 0.289
## F-statistic: 3.439 on 2 and 10 DF, p-value: 0.07302 legend("topleft", p.trans(anova(res)$P[1]), bty = "n")
mtext("Dry versus not dry lake", side = 1, line = 2, outer = TRUE)
Time structured analyses
How does recruitment in lake 2 affect the whole population?
Check how the number of previous life stages correlates with numbers of later life stages.
# Functions for drawing
plot.regr.FUN <- function(x,y, alpha = 0.05, pch = 16, yaxt = "s", lab = names(x)) {
plot(y ~ x, las = 1, pch = pch, yaxt = yaxt)
text(x, y, lab, pos = 3, col = "darkgrey")
res <- lm(y ~ x)
if (anova(res)$"Pr(>F)"[1] <= alpha) {abline(res); legend("topright", p.trans(anova(res)$"Pr(>F)"[1]), bty = "n")}
if (anova(res)$"Pr(>F)"[1] > alpha) {abline(h = mean(y, na.rm = TRUE), lty = 2, col = "darkgrey");legend("topright", "p > 0.05 [NS]", bty = "n")}
}
time.FUN <- function(x,y, time = 5, ytitle = "") {
plot.regr.FUN(x,y)
mtext(ytitle, 2, 3, F, cex = 0.8)
n <- length(x)
for (i in 1:(time-1)) plot.regr.FUN(y = y[-c(1:i)], x = x[-c((n-(i-1)):n)], yaxt = "n")
}
# Temporal autocorrelation in lake 2
data3 <- subset(data2, site == 2 & year <= 2009)
x1 <- data3$eggs; names(x1) <- 1992:2009
x2 <- data3$larvae; names(x2) <- 1992:2009
x3 <- data3$juveniles; names(x3) <- 1992:2009
x4 <- data3$frogs; names(x4) <- 1992:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x1, data3$eggs, ytitle = "Egg masses")## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliabletime.FUN(x2, data3$larvae, ytitle = "Larvae")## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliabletime.FUN(x3, data3$juveniles, ytitle = "Juveniles")## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliabletime.FUN(x4, data3$frogs, ytitle = "Adults and subadults")## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliablemtext("Population in previous years", 1, 1, T, cex = 0.8)
mtext("Temporal autocorrelation in lake 2", 3, 1, TRUE, font = 2)
# Relationships between metrics
pairs(data3[,3:9])
pairs(data2[,10:14])
# Previous life stages effects in lake 2 on the whole bassin
data3 <- subset(data2, site == 2 & year <= 2009)
x1 <- tapply(data2$eggs, data2$year, mean, na.rm = TRUE)[4:21]
x2 <- tapply(data2$larvae, data2$year, mean, na.rm = TRUE)[4:21]
x3 <- tapply(data2$juveniles, data2$year, mean, na.rm = TRUE)[4:21]
x4 <- tapply(data2$frogs, data2$year, mean, na.rm = TRUE)[4:21]
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$eggs, x1, ytitle = "Egg masses")
time.FUN(data3$eggs, x2, ytitle = "Larvae")
time.FUN(data3$eggs, x3, ytitle = "Juveniles")
time.FUN(data3$eggs, x4, ytitle = "Adults and subadults")
mtext("Egg masses in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of eggs masses from lake 2 on the whole basin", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$larvae, x1, ytitle = "Egg masses")
time.FUN(data3$larvae, x2, ytitle = "Larvae")
time.FUN(data3$larvae, x3, ytitle = "Juveniles")
time.FUN(data3$larvae, x4, ytitle = "Adults and subadults")
mtext("Larvae in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of larvae from lake 2 on the whole basin", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$juveniles, x1, ytitle = "Egg masses")
time.FUN(data3$juveniles, x2, ytitle = "Larvae")
time.FUN(data3$juveniles, x3, ytitle = "Juveniles")
time.FUN(data3$juveniles, x4, ytitle = "Adults and subadults")
mtext("Juveniles in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of juveniles from lake 2 on the whole basin", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$frogs, x1, ytitle = "Egg masses")
time.FUN(data3$frogs, x2, ytitle = "Larvae")
time.FUN(data3$frogs, x3, ytitle = "Juveniles")
time.FUN(data3$frogs, x4, ytitle = "Adults and subadults")
mtext("Adults and subadults in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of adults and subadults from lake 2 on the whole basin", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$juveniles[-8], x1[-8], ytitle = "Egg masses")
time.FUN(data3$juveniles[-8], x2[-8], ytitle = "Larvae")
time.FUN(data3$juveniles[-8], x3[-8], ytitle = "Juveniles")
time.FUN(data3$juveniles[-8], x4[-8], ytitle = "Adults and subadults")
mtext("Juveniles in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of juveniles from lake 2 on the whole basin (without 1999)", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$frogs[-8], x1[-8], ytitle = "Egg masses")
time.FUN(data3$frogs[-8], x2[-8], ytitle = "Larvae")
time.FUN(data3$frogs[-8], x3[-8], ytitle = "Juveniles")
time.FUN(data3$frogs[-8], x4[-8], ytitle = "Adults and subadults")
mtext("Adults and subadults in lake 2 in previous years", 1, 1, T, cex = 0.8)
mtext("Effects of adults and subadults from lake 2 on the whole basin (without 1999)", 3, 1, TRUE, font = 2)
The same but only in lake 2:
data3 <- subset(data2, site == 2 & year <= 2009)
x1 <- data3$eggs; names(x1) <- 1992:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x1, data3$eggs, ytitle = "Egg masses")## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliabletime.FUN(x1, data3$larvae, ytitle = "Larvae")
time.FUN(x1, data3$juveniles, ytitle = "Juveniles")
time.FUN(x1, data3$frogs, ytitle = "Adults and subadults")
mtext("Population in previous years", 1, 1, T, cex = 0.8)
mtext("Number of egg masses in previous years", 3, 1, TRUE, font = 2)
x1 <- data3$larvae; names(x1) <- 1992:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x1, data3$eggs, ytitle = "Egg masses")
time.FUN(x1, data3$larvae, ytitle = "Larvae")## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliabletime.FUN(x1, data3$juveniles, ytitle = "Juveniles")
time.FUN(x1, data3$frogs, ytitle = "Adults and subadults")
mtext("Population in previous years", 1, 1, T, cex = 0.8)
mtext("Number of larvae in previous years", 3, 1, TRUE, font = 2)
x1 <- data3$juveniles; names(x1) <- 1992:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x1, data3$eggs, ytitle = "Egg masses")
time.FUN(x1, data3$larvae, ytitle = "Larvae")
time.FUN(x1, data3$juveniles, ytitle = "Juveniles")## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliabletime.FUN(x1, data3$frogs, ytitle = "Adults and subadults")
mtext("Population in previous years", 1, 1, T, cex = 0.8)
mtext("Number of juveniles in previous years", 3, 1, TRUE, font = 2)
x1 <- data3$frogs; names(x1) <- 1992:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x1, data3$eggs, ytitle = "Egg masses")
time.FUN(x1, data3$larvae, ytitle = "Larvae")
time.FUN(x1, data3$juveniles, ytitle = "Juveniles")
time.FUN(x1, data3$frogs, ytitle = "Adults and subadults")## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliable
## Warning in anova.lm(res): ANOVA F-tests on an essentially perfect fit are
## unreliablemtext("Population in previous years", 1, 1, T, cex = 0.8)
mtext("Number of adults and subadults in previous years", 3, 1, TRUE, font = 2)
Lots of plots to figure out whether environmental data explains how many frogs were counted each summer:
# Proportion of occupied lakes as functions of seasonal data
YGF2 <- cbind(matrix(NA, 4,4), YGF[,1:7], rep(NA, 4), YGF[,8:12])
colnames(YGF2) <- 1993:2009
par(mfrow = c(1,1), las = 2, mar = c(4,4,1,1))
barplot(YGF2[,-c(1:4)], beside = TRUE, legend = TRUE, ylim = c(0, 1), args.legend = list(x = "topright", bty = "n"), col = terrain.colors(4), ylab = "Proportion of occupied lakes")
x <- data3$winter[-1]; names(x) <- colnames(YGF2) <- 1993:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, YGF2[1,], ytitle = "Egg masses")
time.FUN(x, YGF2[2,], ytitle = "Larvae")
time.FUN(x, YGF2[3,], ytitle = "Juveniles")
time.FUN(x, YGF2[4,], ytitle = "Adults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Proportion of occupied lakes", 3, 1, TRUE, font = 2)
x <- data3$wc.mean[-1]; names(x) <- colnames(YGF2) <- 1993:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, YGF2[1,], ytitle = "Egg masses")
time.FUN(x, YGF2[2,], ytitle = "Larvae")
time.FUN(x, YGF2[3,], ytitle = "Juveniles")
time.FUN(x, YGF2[4,], ytitle = "Adults")
mtext("Average water content", 1, 1, T, cex = 0.8)
mtext("Proportion of occupied lakes", 3, 1, TRUE, font = 2)
x <- data3$sum.years[-1]; names(x) <- colnames(YGF2) <- 1993:2009
#quartz("", 12, 7)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, YGF2[1,], ytitle = "Egg masses")
time.FUN(x, YGF2[2,], ytitle = "Larvae")
time.FUN(x, YGF2[3,], ytitle = "Juveniles")
time.FUN(x, YGF2[4,], ytitle = "Adults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Proportion of occupied lakes", 3, 1, TRUE, font = 2)
# Winter length, water content (wc) and summer length in lake 2
data3 <- subset(data2, site == 2 & year <= 2009)
x <- data3$winter; names(x) <- 1992:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 2", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 2", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 2", 3, 1, TRUE, font = 2)
# Winter length, water content (wc) and summer length in lake 5
data3 <- subset(data2, site == 5 & year <= 2009)
x <- data3$winter; names(x) <- 1992:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 5", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 5", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 5", 3, 1, TRUE, font = 2)
# Winter length, water content (wc) and summer length in lake 7
data3 <- subset(data2, site == 7 & year <= 2009)
x <- data3$winter; names(x) <- 1992:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 7", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 7", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 7", 3, 1, TRUE, font = 2)
# Winter length, water content (wc) and summer length in lake 20
data3 <- subset(data2, site == 20 & year <= 2009)
x <- data3$winter; names(x) <- 1992:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 20", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 20", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 20", 3, 1, TRUE, font = 2)
# Winter length, water content (wc) and summer length in lake 4
data3 <- subset(data2, site == 4 & year <= 2009)
x <- data3$winter; names(x) <- 1992:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 4", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 4", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 4", 3, 1, TRUE, font = 2)
# Winter length, water content (wc) and summer length in lake 1
data3 <- subset(data2, site == 1 & year <= 2009)
x <- data3$winter; names(x) <- 1992:2009
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(x, data3$eggs, ytitle = "Egg masses")
time.FUN(x, data3$larvae, ytitle = "Larvae")
time.FUN(x, data3$juveniles, ytitle = "Juveniles")
time.FUN(x, data3$frogs, ytitle = "Adults and subadults")
mtext("Winter season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 1", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$wc.mean, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$wc.mean, data3$larvae, ytitle = "Larvae")
time.FUN(data3$wc.mean, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$wc.mean, data3$frogs, ytitle = "Adults and subadults")
mtext("Water content (yearly average)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 1", 3, 1, TRUE, font = 2)
par(mfrow = c(4,5), mar = c(2,0,0,0), oma = c(2,4,3,1))
time.FUN(data3$sum.years, data3$eggs, ytitle = "Egg masses")
time.FUN(data3$sum.years, data3$larvae, ytitle = "Larvae")
time.FUN(data3$sum.years, data3$juveniles, ytitle = "Juveniles")
time.FUN(data3$sum.years, data3$frogs, ytitle = "Adults and subadults")
mtext("Summer season length (relative to year)", 1, 1, T, cex = 0.8)
mtext("Life stages in lake 1", 3, 1, TRUE, font = 2)
Population size vs. lake dried out or not
table(water$Year)##
## 2002 2003 2005 2006
## 49 29 41 30drought.prob <- data.frame(min_surf = as.numeric(min_surf), year = as.numeric(rep(rownames(min_surf), times = dim(min_surf)[2])), site = as.numeric(rep(colnames(min_surf), each = dim(min_surf)[1])))
#--------------------------------------------------------
# Lake dried (0), or not (1)
#--------------------------------------------------------
drought.prob$drought <- drought.prob$min_surf
drought.prob$drought[drought.prob$drought > 0] <- 1
drought.prob2 <- drought.prob
drought.prob2$year <- drought.prob2$year+1
#--------------------------------------------------------
# Match population summaries (data2) with drought prob
# The year during and after the drought
#--------------------------------------------------------
drought.prob <- merge(data2, drought.prob)
drought.prob2 <- merge(data2, drought.prob2)
#--------------------------------------------------------
# Compare populations between lakes that dried versus not dried
#--------------------------------------------------------
table(drought.prob$drought, drought.prob$site, drought.prob$year)## , , = 2002
##
##
## 1 2 3 4 5 6 7 8 10 20 22
## 0 0 1 0 0 0 1 0 1 1 0 0
## 1 1 0 1 1 1 0 1 0 0 1 0
##
## , , = 2003
##
##
## 1 2 3 4 5 6 7 8 10 20 22
## 0 0 0 0 0 0 1 0 1 1 0 0
## 1 1 1 1 1 1 0 1 0 0 0 1
##
## , , = 2005
##
##
## 1 2 3 4 5 6 7 8 10 20 22
## 0 0 0 0 0 0 0 0 1 1 0 0
## 1 1 1 1 1 1 1 1 0 0 1 1
##
## , , = 2006
##
##
## 1 2 3 4 5 6 7 8 10 20 22
## 0 0 1 0 0 0 0 0 0 0 0 0
## 1 1 0 1 1 1 1 1 1 1 1 0par(mfrow = c(2,4), mar = c(1,2,2,1), oma = c(4,3,1,0))
boxplot(I(eggs+1) ~ drought, drought.prob, log = "y", main = "Egg masses", las = 1)
res <- lmer(log(eggs+1) ~ drought + (1|site), drought.prob); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(eggs + 1) ~ drought + (1 | site)
## Data: drought.prob
##
## REML criterion at convergence: 94
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4924 -0.3973 -0.2604 -0.1466 2.2097
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.1691 0.4113
## Residual 0.4563 0.6755
## Number of obs: 41, groups: site, 11
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.40387 0.27979 24.47525 1.444 0.162
## drought 0.06342 0.29716 37.59510 0.213 0.832
##
## Correlation of Fixed Effects:
## (Intr)
## drought -0.811 legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Pop. size (during drought)", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ drought, drought.prob, log = "y", main = "Larvae", las = 1)
res <- lmer(log(larvae +1) ~ drought + (1|site), drought.prob); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(larvae + 1) ~ drought + (1 | site)
## Data: drought.prob
##
## REML criterion at convergence: 147.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.87044 -0.44350 -0.08483 0.48452 2.55453
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 6.919 2.6303
## Residual 0.959 0.9793
## Number of obs: 41, groups: site, 11
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.3667 0.9001 14.5264 4.851 0.000231 ***
## drought -1.1922 0.5132 31.1300 -2.323 0.026879 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## drought -0.440 legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(juveniles+1) ~ drought, drought.prob, log = "y", main = "Juveniles", las = 1)
res <- lmer(log(juveniles +1) ~ drought + (1|site), drought.prob); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(juveniles + 1) ~ drought + (1 | site)
## Data: drought.prob
##
## REML criterion at convergence: 142.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.60978 -0.28228 0.03419 0.55254 1.49806
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 2.574 1.604
## Residual 1.130 1.063
## Number of obs: 41, groups: site, 11
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.4202 0.6587 20.6056 2.156 0.0431 *
## drought 0.8475 0.5373 34.7451 1.577 0.1238
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## drought -0.628 legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(frogs+1) ~ drought, drought.prob, log = "y", main = "Adults and subadults", las = 1)
res <- lmer(log(frogs +1) ~ drought + (1|site), drought.prob); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(frogs + 1) ~ drought + (1 | site)
## Data: drought.prob
##
## REML criterion at convergence: 88
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.56119 -0.62143 0.05131 0.64305 1.51461
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 1.2259 1.1072
## Residual 0.2258 0.4751
## Number of obs: 41, groups: site, 11
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.1267 0.3920 15.8250 5.426 5.82e-05 ***
## drought 0.4606 0.2475 31.7976 1.861 0.072 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## drought -0.487 legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(eggs+1) ~ drought, drought.prob2, log = "y", main = "", las = 1)
res <- lmer(log(eggs+1) ~ drought + (1|site), drought.prob2); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(eggs + 1) ~ drought + (1 | site)
## Data: drought.prob2
##
## REML criterion at convergence: 73.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3109 -0.6362 -0.3160 0.5727 1.9577
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.171 0.4135
## Residual 0.782 0.8843
## Number of obs: 27, groups: site, 11
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.2778 0.3920 17.3623 0.709 0.488
## drought 0.4815 0.4519 19.1957 1.066 0.300
##
## Correlation of Fixed Effects:
## (Intr)
## drought -0.840 legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Pop. size (after drought)", side = 2, line = 3, outer = FALSE)
boxplot(I(larvae+1) ~ drought, drought.prob2, log = "y", main = "", las = 1)
res <- lmer(log(larvae +1) ~ drought + (1|site), drought.prob2); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(larvae + 1) ~ drought + (1 | site)
## Data: drought.prob2
##
## REML criterion at convergence: 115.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.08430 -0.44864 -0.00434 0.34852 1.95495
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 2.770 1.664
## Residual 3.275 1.810
## Number of obs: 27, groups: site, 11
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.185 1.024 17.991 2.135 0.0468 *
## drought 1.785 1.137 22.835 1.570 0.1303
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## drought -0.798 legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(juveniles+1) ~ drought, drought.prob2, log = "y", main = "", las = 1)
res <- lmer(log(juveniles +1) ~ drought + (1|site), drought.prob2); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(juveniles + 1) ~ drought + (1 | site)
## Data: drought.prob2
##
## REML criterion at convergence: 99.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0094 -0.4241 -0.0583 0.4846 1.8962
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 2.006 1.416
## Residual 1.536 1.239
## Number of obs: 27, groups: site, 11
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.5573 0.7783 19.2347 0.716 0.483
## drought 1.9496 0.8404 24.4645 2.320 0.029 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## drought -0.773 legend("topleft", p.trans(anova(res)$P), bty = "n")
boxplot(I(frogs+1) ~ drought, drought.prob2, log = "y", main = "", las = 1)
res <- lmer(log(frogs +1) ~ drought + (1|site), drought.prob2); summary(res)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(frogs + 1) ~ drought + (1 | site)
## Data: drought.prob2
##
## REML criterion at convergence: 75.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.97522 -0.41394 -0.00719 0.24147 2.13657
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.5121 0.7156
## Residual 0.6661 0.8161
## Number of obs: 27, groups: site, 11
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.3769 0.4518 18.5096 3.048 0.00677 **
## drought 1.7862 0.5043 22.7391 3.542 0.00176 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## drought -0.802 legend("topleft", p.trans(anova(res)$P), bty = "n")
mtext("Dry versus not dry lake", side = 1, line = 2, outer = TRUE)